Claudio Iturra, 2023

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from iwaves.kdv.vkdv import  vKdV as KdV
In [2]:
plt.rcParams['axes.labelsize'] = 'medium'
plt.rcParams['font.size'] = 14
In [3]:
def depth_tanh(beta, x):
    """
    Hyperbolic tangent shelf break

    H - total depth
    h0 - shelf height
    x0 - shelf break x location
    lt - shelf break width
    """
    
    H, h0, x0, lt = beta

    return H-0.5*h0*(1+np.tanh((x-x0)/(0.5*lt)))
In [6]:
# Model parameters
N = 3000 # Number of grid points
dx = 50. # horizontal grid spacing

H = 600
h0 = 450
bathy_params = [H, h0, 75e3, 70e3]  # See the depth_tanh function for the mean

# The horizontal domain
x = np.arange(0, N*dx, dx)

# The depth
h = depth_tanh(bathy_params, x) 

plt.figure()
plt.plot(x, -h)
plt.title('KdV model bathymetry')
plt.ylabel('Depth [m]')
plt.xlabel('x [m]')
plt.ylim(-H, 0)
Out[6]:
(-600.0, 0.0)
In [7]:
def rho_double_tanh(beta, z):
    """
    Double hyperbolic tangent density profile model
    """
    return beta[0] - beta[1]*(np.tanh((z+beta[2])/beta[3])
        + np.tanh((z+beta[4])/beta[5]))
In [8]:
# See the references above for the mean of the
rho_params =[1023.68, # Mean density
     1.22, # Density jump
     156.7, # Depth of the first pycnocline
     53.6,# Width of the first pycnocline
     73.1,# Depth of the first pycnocline
     40.2, # Width of the second pycnocline
        ] 

# Number of vertical levels
Nz = 50

z = np.linspace(-H,0,Nz)
rhoz = rho_double_tanh(rho_params,z)

plt.figure()
plt.plot(rhoz, z)
plt.ylabel('Depth [m]')
plt.xlabel('Density [kg m$^{-3}$]')
Out[8]:
Text(0.5, 0, 'Density [kg m$^{-3}$]')
In [9]:
# IMEX options (These are weights for the numerical time integration scheme)
imex={
        'MCN_AX2':(1/8., 3/8.),
        'AM2_AX2':(1/2., 1/2.),
        'AI2_AB3':(3/2., 5/6.),
        'BDF2_BX2':(0.,0.),
        'BDF2_BX2s':(0.,1/2.),
        'BI2_BC3':(1/3.,2/3.),
}
In [10]:
dt = 15.
mode = 0 # Mode=0 corresponds to mode-1 waves
imexscheme = 'AM2_AX2'
c_im = imex[imexscheme][0]
b_ex = imex[imexscheme][1]

kdvargs = dict(
   N=N,
   dx=dx,
   dt=dt,
   spongedist = 5e3,
   spongetime = 360.,
   Nsubset = 10,
   nonhydrostatic=1.,
   nonlinear=1.,
   c_im=c_im,
   b_ex=b_ex,
   verbose=True,
)

## Initialise the class
mykdv = KdV(rhoz, z, h, x, mode, **kdvargs)

Calculating eigenfunctions...
0.0 % complete...
5.0 % complete...
10.0 % complete...
15.0 % complete...
20.0 % complete...
25.0 % complete...
30.0 % complete...
35.0 % complete...
40.0 % complete...
45.0 % complete...
50.0 % complete...
55.0 % complete...
60.0 % complete...
65.0 % complete...
70.0 % complete...
75.0 % complete...
80.0 % complete...
85.0 % complete...
90.0 % complete...
95.0 % complete...
In [11]:
# Plot some of the environmental parameters as a function of distance
plt.figure(figsize=(6,10))
ax1=plt.subplot(411)
ax1.plot(mykdv.x, mykdv.c)
plt.ylabel('c [m s$^{-1}$]')
ax1.set_xticklabels([])
ax2=plt.subplot(412)
ax2.plot(mykdv.x, mykdv.alpha)
plt.ylabel(r'$\alpha$ [s$^{-1}$]')
ax2.set_xticklabels([])

ax3=plt.subplot(413)
ax3.plot(mykdv.x, mykdv.beta)
plt.ylabel(r'$\beta$ [m$^3$ s$^{-1}$]')
ax3.set_xticklabels([])

ax2=plt.subplot(414)
ax2.plot(mykdv.x, 1/mykdv.Q)
plt.ylabel(r'Linear amplification "Q"')
Out[11]:
Text(0, 0.5, 'Linear amplification "Q"')
In [12]:
# Plot the density profile and vertical eigenfunction
plt.figure(figsize=(6,6))
ax1=plt.subplot(121)
ax1.plot(mykdv.Phi[:,0],mykdv.Z[:,0], )
plt.xlabel('$\phi$')

ax2=plt.subplot(122)
ax2.plot(mykdv.rhoZ[:,0],mykdv.Z[:,0], )
plt.xlabel('Density [kg m$^{-3}$]')
ax2.set_yticklabels([])
ax1.set_ylabel('Depth [m]')
Out[12]:
Text(0, 0.5, 'Depth [m]')
In [13]:
# Plot the density profile and vertical eigenfunction (upstream)
plt.figure(figsize=(6,6))
ax1=plt.subplot(121)
ax1.plot(mykdv.Phi[:,2000],mykdv.Z[:,2000], )
plt.xlabel('$\phi$')

ax2=plt.subplot(122)
ax2.plot(mykdv.rhoZ[:,2000],mykdv.Z[:,2000], )
plt.xlabel('Density [kg m$^{-3}$]')
ax2.set_yticklabels([])
ax1.set_ylabel('Depth [m]')
Out[13]:
Text(0, 0.5, 'Depth [m]')
In [14]:
def bcfunc(a0, period, t):
    omega = 2*np.pi/period
    return a0*np.sin(omega*t)
In [20]:
# Specify the runtime in seconds
runtime = 2*86400 # 3 days
ntout = 600//dt # how often to output the model solution

# Boundary condition parameters
a0 = 25.
period = 20*3600.

# Number of time steps to run
nsteps = int(runtime//dt)

# Plot up the boundary function
t = np.arange(0, nsteps*dt, dt)
plt.figure()

plt.plot(t, bcfunc(a0, period, t))
plt.ylabel('A [m]')
plt.xlabel('t [s]')
Out[20]:
Text(0.5, 0, 't [s]')
In [23]:
t=[]
A=[]
mykdv.t=0
for ii in range(nsteps):
    if ii%ntout==0:
        # Store snapshots of the solution every "ntout" steps
        t.append(mykdv.t)
        A.append(mykdv.B_n_p1*1) # the solution is stored as the "B_n_p1" 
        
    if mykdv.solve_step(bc_left=bcfunc(a0,period,mykdv.t)) != 0:
        print('Blowing up at step: %d'%ii)
        break
        

        
t = np.array(t)
A = np.array(A)
In [24]:
#"x-t"/Hovmoller plot of the solution

plt.figure()
plt.pcolormesh(x, t, A, cmap='RdBu')
plt.ylabel('t [s]')
plt.xlabel('x [m]')
Out[24]:
Text(0.5, 0, 'x [m]')
In [25]:
# Animate the solution
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
In [26]:
tstep = 0

fig=plt.figure(figsize=(12,5))
ax=plt.subplot(111)
ll, = ax.plot(x, A[tstep,:])
plt.xlim(x[0],x[-1])
plt.ylabel('A [m]')
plt.xlabel('x [m]')
txt= plt.text(0.05, 0.9, "t={} [s]".format(tstep*mykdv.dt*ntout), transform=ax.transAxes)
plt.ylim(-3*a0, 3*a0)

def update(tstep):

        ll.set_ydata(A[tstep,:])
        txt.set_text("t={} [s]".format(tstep*mykdv.dt*ntout))
        
        return txt
    
anim = FuncAnimation(fig, update, frames=range(0, t.shape[0]))
HTML(anim.to_jshtml())
Out[26]:
In [ ]: